Slow quench dynamics of Mott-insulating regions in a trapped Bose-gas 
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We investigate the dynamics of Mott-insulating regions of a trapped bosonic gas as the interaction strength 
is changed linearly with time. The bosonic gas considered is loaded into an optical lattice and confined to a 
parabolic trapping potential. Two situations are addressed: the formation of Mott domains in a superfluid gas 
as the interaction is increased, and their melting as the interaction strength is lowered. In the first case, de- 
pending on the local filling, Mott-insulating barriers can develop and hinder the density and energy transport 
throughout the system. In the second case, the density and local energy adjust rapidly whereas long range 
correlations require longer time to settle. For both cases, we consider the time evolution of various observ- 
ables: the local density and energy, and their respective currents, the local compressibility, the local excess 
energy, the heat and single particle correlators. The evolution of these observables is obtained using the time- 
dependent density-matrix renormalization group technique and comparisons with time-evolutions done within 
the Gutzwiller approximation are provided. 

PACS numbers: 05.70.Ln, 02.70.-c, 05.30.Rt, 67.85.Hj 



I. INTRODUCTION 

Due to their good isolation from the environment and to 
their tunability, ultra-cold quantum gases are ideal candidates 
to explore systems away from equilibrium 1 1 ] . Cold atoms are 
well suited to explore situations where the Hamiltonian de- 
scribing a system is slowly varied with time. Understanding 
the physical implications of such slow quenches is of great 
theoretical and practical importance to shed light on the co- 
herent evolution of quantum systems and to devise methods 
to prepare complex quantum phases. Seminal works on the 
dynamics of classical systems near a second order phase tran- 
sition conducted by Kibble [2] and Zurek |3] identified that 
the defect production rate as a function of the ramp velocity is 
described by a scaling law when the system crosses a critical 
point. However, despite many recent theoretical advances |4- 
13111 . the response of strongly-correlated quantum gases to the 
slow quench of a Hamiltonian parameter is still far from be- 
ing fully understood. Meanwhile, on the experimental side, 
considerable efforts have been devoted to understand the dy- 
namics of interacting bosonic atoms when the depth of the 
optical lattice is varied ll32l - [36ll or when the slow quench of 
an effective parameter is performed IstI [ssll . 

In relation to these experimental protocols, in Ref. [13, HH 
[39I and 113, the presence of a parabolic trapping potential 
was found to significantly influence the dynamics. Two dy- 
namical regimes have been shown to exist when interacting 
atoms loaded into an optical lattice and confined to a trap are 
subjected to a slow change of the interaction strength. For 
short ramp times, the evolution is dominated by intrinsic lo- 
cal dynamics, which is also present in a homogeneous system, 
whereas, for longer ramp times the density redistribution can 
play an important role. 

In this work, we study the response of bosonic atoms to 
a linear change of the interaction strength. As these atoms 



are confined to one-dimensional tubes and loaded into an op- 
tical lattice running along the tubes main axis, the physics 
for a wide range of parameters is well described by the one- 
dimensional Bose-Hubbard model. Here our main objective is 
to understand the evolution, as a function of the ramp time, of 
the local and non-local observables of the quantum gas. We 
focus on the crucial issues of the formation and melting of 
Mott domains and on how the adiabatic limit is approached. 



The article is structured as follows: In Sec. Ill Al we in- 
troduce the model and the time-dependent protocol. Sec- 
tions HlB] and HTC] detail the methods and the theoretical def- 
initions. These two sections can be skipped by readers more 
interested in the main phenomena. In Sec. |III1 we turn to the 
description of the evolution resulting from the increase of the 
interaction strength. In Sec. IIII A 1 1 we focus on the occur- 
rence of two dynamical regimes, the intrinsic dynamics and 
the dynamics induced by the trapping potential, and explain 
their origin (Sec. IIII A 21 ). Afterwards, we direct our atten- 
tion to the formation of "Mott barriers" which strongly block 
the equilibration process (Sec. IIII A 31 ). the energy transport 
(Sec. Ill C 4l i and the evolution of longer range correlations 
(Sec. IIII Bt . Then, in Sec. |IV| we consider the opposite case 
of melting the Mott domains occurring when the interaction 
strength is lowered and, in particular, we point out the long 
equilibration times for long range correlation functions. For 
both situations, we characterize the time evolution consider- 
ing various observables such as the density, the compressibil- 
ity, the energy, various particle correlators, and the momen- 
tum distribution which is related to the interference patterns in 
time-of-flight measurements. These results are supplemented 
by a detailed analysis of the physical mechanisms responsible 
for the presence of the intrinsic and global dynamical regimes. 
We further show that the different time-scales can be identi- 
fied experimentally from interference patterns. Our numerical 
results are obtained from the time-dependent density-matrix 
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renormalization group method (t-DMRG). We also compare 
these quasi-exact resuhs to time-evolutions done within the 
mean-field Gutzwiller method (Sec. IIII A 5b to identify the 
limitations of the latter approach and pinpoint the qualitative 
physical insights it provides. The present work extends sub- 
stantially our previous results on the same setup iITtIi . 

II. MODEL AND CONSERVATION LAWS 

A. Hamiltonian and time-dependent protocol 

Ultracold bosons in optical lattices are, in a wide parameter 
r^ime, well described by the Bose-Hubbard Hamiltonian ll4ll 

/ / I 

with the operator creating a boson at site / and fii — bjbi 
the local density operator. The total number of atoms is fixed 
to N. The first term of the Hamiltonian corresponds to the 
kinetic energy of the atoms with the hopping amplitude J 
and the second to the onsite interaction of strength U. The 
site-dependent chemical potential /i; accounts for an external 
confinement. We consider here a one-dimensional geometry 
(tube) and use either a homogeneous (ni = 0) or a harmonic 
trapping potential of the form /i; = —Vo{l — {L + l)/2)^, 
with L the number of sites in the tube (open-boundary condi- 
tions) used in our simulations. We assume an experimentally 
realistic strength for the trapping potential of Vb = 0.006J 
and particle number N = 24, 48. For these parameters 
the choice L = 64 assures that edge effects are not impor- 
tant. One non-trivial aspect of the model is that it is non- 
integrable ll43il44ll for non-zero J and U. Further, at commen- 
surate fillings, a quantum phase transition from a superfluid 
to a Mott-insulating state occurs (at {U/J)c ~ 3.4 for unity 
filling in one-dimension |45, 46]). This phase transition is ac- 
companied by the opening of a gap in the low-energy excita- 
tion spectrum, which strongly modifies the ground-state, ther- 
modynamic and transport properties. At incommensurate fill- 
ings a crossover between a superfluid and a Tonks-Girardeau, 
or hard-core boson, gas occurs in equilibrium. This distinct 
behavior at commensurate and incommensurate fillings im- 
plies that in a trapped system, different states can coexist in 
spatially separated regions |47-49]. For instance, for strong 
enough interaction, a Mott-insulating plateau with commen- 
surate filling, surrounded by a superfluid region, emerges. 

Regarding the time-dependent protocol, we consider a slow 
quench of the interaction strength U{t) which can be achieved 
experimentally using a suitable Feshbach resonance llsoll . Dif- 
ferent time-dependent protocols have been considered in pre- 
vious works in homogeneous systems using several analyti- 
cal or numerical approximation schemes iH |^ [l^ [l^ |20| - 
112, HH]. For sake of simplicity and generality, the variation 
in time is chosen to be linear, starting from Ui up to a fi- 
nal value Uj: U{t) = Ui + ^5U with r the ramp time and 
5U = Uf — Ui the quench amplitude. The real-time evolution 



starts from the ground state corresponding to Ui. The labels 
i/ f are used for the initial and final ground state values, re- 
spectively. The limit t — > 0, i.e. the sudden quench limit, has 
been studied intensively in the Bose-Hubbard model using an- 
alytical |i23,,51,,52J and numerical methods |i53n5^. 



B. Methods 

1. t-DMRG 

Accurate ab-initio numerical simulations of the time evo- 
lution of the quantum gas are carried out using the t-DMRG 
technique |'60'-'64] . The time-evolution is implemented using 
the second order Trotter- Suzuki decomposition. The dimen- 
sion of the effective space is a few hundred states and the time- 
step is adjusted with the ramp velocity. We introduce a cutoff 
value of M = 5 or 6 in the number of onsite bosons as higher 
boson occupancies are negligible in the situations considered 
here. 



2. Gutzwiller variational method 

In this section, we present how to determine the evolu- 
tion of the system within the Gutzwiller mean-field method 
|65, 66]. This approximation has been used before to describe 
the evolution during a slow change of the lattice depth in a 
higher dimensional trapped Bose-Hubbard system ll40l l67ll . 
The Gutzwiller method is based on a variational ansatz of 
the many-body wave-function Y^) = [X^n, 
where \ni) is the Fock state on site I with rii particles and 
c;_„, are the variational parameters. The ground state for a 
given Hamiltonian is obtained by minimizing the total energy 
Egw'- 

EG^^^jY{{bi)*{bi+i)+c.c^ 

l.rii 1,111 

where (S/) = X)„, V^i + 1 cf „^Ci,„,+i and * denotes com- 
plex conjugation. The validity of the Gutzwiller method 
in evaluating static observables of one-dimensional systems 
has been studied, for example, in Ref. The Gutzwiller 
approach predicts, in one-dimension, a phase transition at 
{U/J)c = 2 (1 + V2)2 ~ 11.7 16?]. The superfluid phase 
is signaled by a non-vanishing order parameter (6;) and, for 
a small interaction strength, the properties of local quanti- 
ties are reasonably well approximated. In contrast, the Mott- 
insulating phase is characterized by a vanishing order param- 
eter and vanishing local compressibility, thus neglecting com- 
pletely particle fluctuations which are present in the real Mott- 
insulating phase. 

The time evolution for the coefficients c;^„, (t) can be read- 
ily derived [661 from the Schrodinger equation. The equations 
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are 

ihdtci^mit) = - 1) - fJ-ini^ ci^mit) 

~jVWTT[{bi-i)* + (bi+i)*] Q,„,+i(t) 

-Jy/rTiiibi-i) + (h+i)] ci.ni-i{t) . 

They are solved numerically by implementing a split-step 
method. 



C. Observables and definitions from continuity equations 

1. Correlations and interference pattern 

We define the one-body correlation function between sites 
I and m as 

9i,m = ^{b]br,i + h-c.) . (1) 

While g is not easily accessible, time-of-flight techniques 
measure an interference pattern related to the momentum 
distribution of the correlated gas. Neglecting the Wannier- 
function envelope, the interference pattern is given by 

l,m 

where a is the lattice spacing. For a superfluid state, N{k) 
is expected to be strongly peaked around zero momentum, 
whereas for a Mott-insulating state the interference pattern 
should be rather flat. 



2. General expression for the continuity equation 

In this section, we want to determine, in the Schrodinger 
picture, the current operators corresponding to a given observ- 
able 0{t) = {i}{t)\d{t)\-4}{t)), where d{t) can explicitly de- 
pend on time, using the associated continuity equation. In the 
following, the shorthand notation (• • •) stands for the expec- 
tation value ■ a given time. The continuity 
equation takes, on general grounds, the following form: 

hdtO{t) = -div(J°) + {S") . (3) 

is the current operator for which we have 

= i{[H{t),d{t)]) (4a) 

= -{{J?,i+i)-{J?-u))- (4b) 

The second equality is the specialization to a one-dimensional 
lattice for an observable located around site /, with incom- 
ing and outgoing currents. As we are interested in studying a 
one-dimensional system, ( l4b] ) is used throughout. The source 
operator S'°(i) = hdtd{t) is non-zero only for an explicitly 
time-dependent operator. 



Interestingly, integrating (O between times and r, taking 
the adiabatic limit r — > cx) and doing the change of variable 
t ^ U (in the source term integral) enables one to express 
the integrated contribution of currents only as a function of 
ground state expectation values 

oo Uf 

J ^div(jO)(t) = o,-0; + Jdu {Mu)\dud\Mu)) , 

(7. 

(5) 

where \iIjq{U)) is the ground state corresponding to U. The 
integral on the right-hand side is taken along the adiabatic 
path. This remark is important from a numerical perspec- 
tive because the right-hand side can be efficiently computed 
via ground state techniques while the left-hand side would re- 
quire time-dependent simulations over very long times, which 
is not feasible. 

As explained in the introduction, the main objective of this 
work is to characterize how particles and energy redistribute 
when the interaction strength, U (t), is ramped up or down. 
Therefore, we introduce below the relevant quantities and 
the physically significant terms of their associated continu- 
ity equations. A few commutators, useful in the derivation of 
these continuity equations, are provided in AppendixlAl 

3. Local observables on sites and bonds 

As a first example, ([3]) can be used to derive the particle 
current associated with the local density n/: 

Xk = j;%^iJ{bibi~blbk). (6) 

This current is defined between sites / and k and there is no 
source term associated with ni as it is not explicitly time- 
dependent. As the particle current appears quite often in the 
rest of this article, from now on, it will be denoted as ji^k- 
Finally, it is instructive to note that for a homogeneous and 
translationally invariant system, the local density is constant 
at all time due to the conservation of the total number of par- 
ticles. 

In order to better understand the different time-scales in- 
volved during the evolution, it is also useful to consider sep- 
arately the evolution equation for the density fluctuations nf. 
This quantity is essential to our comprehension of the Bose- 
Hubbard model and is related to the the local compressibility 
K; = (nf) — (fii)'^. Using ^ and after some algebra, we find 
that the evolution of is controlled by a "density-assisted" 
or "correlated" current 

Ji^l+i = fiiin+i + ji,i+ini (7) 

(note that J^l^^ j = niji-i.i + ji-i.ini)- The origin of this 
density-assisted current, mixing j and vP operators, comes 
from the evolution equation for the onsite occupancy proba- 
bilities discussed in Appendix IB] In equilibrium, in our sys- 
tem, the average of (|7]) computed in the ground state vanishes, 
as for the particle current operator 
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The same strategy is used to get the current operators asso- 
ciated with observables living on bonds, such as the local ki- 
netic energy operator Ki^i+i (or nearest-neighbor one-particle 
correlation) defined by 



Ki,k - -Jiblk + h\hk) 



(8) 



between site / and k. We find that the incoming current asso- 
ciated with Ki^i+i reads 



J, 



uit) . 



i+i 



jLl+lfll) 



(9a) 
(9b) 

(9c) 



showing the interplay of the correlated and usual particle cur- 
rents. It is worth noticing that for a homogeneous system, the 
evolution of local kinetic fluctuations is directly related to that 
of the density fluctuations since in this case 



(10) 



Thus, even in the homogeneous limit, the time dependence of 
U{t) affects the evolution of the local kinetic term or nearest- 
neighbor correlations. Eq. ( fTOl i is also straightforwardly ob- 
tained from the evolution of the total energy discussed be- 
low. Note, this equation is not valid for inhomogeneous gases 
where the balance of particle currents can be non-zero. 

Similarly, the current operator associated with the parti- 
cle current itself contains density-assisted hoppings, follow- 
ing the expression 



pi.i+i 



=^J■lK, 



1,1+1 

U(t) 



(11a) 
(lib) 
(11c) 

(lid) 



We also give for clarity the outgoing current opera- 
tor: j/|i7+2 = fM + lKlU + l + IJ^fli + i + jki^i+2 

-^^{ni+iki^i+i +ki^i+ihi+i). It is worth mentioning that 
the correlated current and hopping terms in (|9]l and (fTTT i all 
come with the interaction strength as a prefactor and disap- 
pear for a non-interacting gas. Their behavior is thus strongly 
affected by the presence of interactions. Finally, as in the next 
section the time-derivative of the particle current will be of 
great use to understand the mechanisms responsible for the 
evolution of the density profile, we provide here its full ex- 
pression: 

fidt{ji,i+i) - ^J'l+l){klJ+l) (12a) 
+ 2j2((n,)-(n,+i)) (12b) 
+ Ji{ki.i,i+i) - {ki^i+2)) (12c) 

- ^{{ni ~ ni+i)ki^i+i + h.c.) . (12d) 



4. Energy and heat 

We now turn to the transport of energy by first defining the 
bond-symmetric local energy operator as 



1 



+ + U[t)li - , (13) 



where = — l)/2 is the operator related to the interac- 
tion energy. In this case, we find that the energy current J^l^ ; 
is given by 



l-l.l 



-Jl-l,l 



J 



{jl-2,l + jl-lj+l) 



2 

Uit) 



(14a) 
(14b) 
(14c) 



(n;-! +ni)ji^ij +ji-i,iini-i + m) , 

(14d) 



in which we naturally recover the particle and correlated cur- 
rents appearing in (|6]l, dTji and (|9]l. In addition, since the en- 
ergy operator is explicitly time-dependent and therefore not a 
conserved quantity during the protocol, we have the following 



source term 



hdtUit)ii 



(15) 



which shows the importance of the density fluctuations in the 
energy production. In particular, the total energy £'(i) = 
(Hit)) satisfies the relation 

dtEit) = {midtUmt)) - [dtUit)]Y,{Ii){t) , (16) 



i.e., the energy put in the system is directly related to the evo- 
lution of the density fluctuations. For an inhomogeneous sys- 
tem, there are two contributions to the local energy production 
as seen from ( fT4l i and (fTsT i: one from currents and correlated 
currents and one from the external driving of the system. Sum- 
ming up the total energy, the contribution from currents must 
vanish to fulfill (fTSI l. but locally, one may have energy redis- 
tribution. We can define the heat produced in the system as 
the energy of the atoms at the final time compared to that of 
the ground state for the final interaction strength 



Q(t) - Eir) - 



SU 

Eo,.— Eoj + — dtV(/;)W, (17) 



with Ef) i/f the ground state energies. Note that (/;) is ac- 
cessible experimentally, which makes it possible to measure 
the interesting Qir) dependence. We can quickly check 
that this formula gives back the correct results in the sudden 
quench and adiabatic quench limits. In the sudden quench 
limit, \i)it)) = \iJoiUi)) which yields Q(0) = Eq^, - 
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Eqj + SUJ2i{^i)o,i- This means that the heat only de- 
pends on ground state properties of the corresponding ini- 
tial and final parameters. In the adiabatic case, we have 
\ijj{t)) — \ipo{U{t))) along the adiabatic path so the integral 

can be reexpressed as Jjj^ dU '^i{Ii)vi{U), with {Ii)o{U) = 

{ipo{U)\Ii\ipo{U)) . Using Feynman-Hellmann theorem over 
U, it is clear that this integral cancels Eo^i — Eqj to make 
Q(oo) = 0. 

One can define a local excess energy qi as the difference in 
local energies between the final energies and the ground state 
expectation for the final parameters: 

qiiT)^{hi)iT)~{hi)oj . (18) 

The local excess energy produced splits up into three different 
contributions 

qiiT)^{hi)o,r~{hi)oj (19a) 

- - / dtdw{j''){t) (19b) 
" Jo 

+— r dt{i,){t), (19c) 

T Jo 

where the first term is simply the local ground state energies 
difference (independent of r), the second term is the inte- 
grated contribution of energy currents, and the last term is the 
integrated contribution due to the external operator While 
Q{t) is necessarily non-negative, <7;(t) can be negative or 
positive depending on the relative contributions of each term. 

Finally, using Q with O ^ hi allows one to calculate these 
quantities in the adiabatic limit; 

oo Uf 

J jdiw{j'''){t)^{h,)o^;-{h,)oj+ jdU {ii)o{U), 

Ui 

(20) 

where the right-hand side can be computed accurately using 
numerical techniques. 

With this set of equations in mind, we are now ready to 
identify the different driving forces responsible for the system 
evolution when the interaction strength is ramped up or down. 

III. DIGGING A MOTT DOMAIN IN A SUPERFLUID 

A. Evolution of local quantities from t-DMRG 

In this section, we consider a linear quench from Ui — 4 J 
to Uf = Q J. Ui is close to the homogeneous superfluid- 
Mott transition point and Uf lies deeper in the Mott-insulating 
regime. We compare two typical situations: (i) the number of 
particles is chosen low enough in order for the maximal filling 
to remain below unity at all times (N — 24); (ii) N is suffi- 
ciently large so that, at Uf, the corresponding ground state 
density profile has a Mott-insulating "shell" and a superfluid 
center (N = 48). We focus on different aspects of the dynam- 
ics: time-scales, role of insulating domains on particle trans- 
port, energy production and transport, and their experimental 
signature. 



N = 24 N = 48 



at site I = 32 at site I = 32 




X [h/J] T [h/J] 



FIG. 1. (color online). Slow quench from Ui — lo Uf = 6 J. 
Evolution of local observables in the presence of a trap as a func- 
tion of the ramp time r, and compared with that of a homogeneous 
system (open symbols) having the same initial local density. Observ- 
ables are the density ni, compressibility k;, occupancy probabilities 
Po and Pi, neighboring correlation gi^i+i and the particle current 
ji,t+i = Subplots correspond to two different total number 

of particles N — 24 and 48, and two different sites: I — 18 and the 
central site 1 — 32 (cf. Fig.|2]for the location of these sites). 



1. Existence of two dynamical regimes 

In Fig. [T] the final values (t = t) of most of the local ob- 
servables introduced before (density, local compressibility, lo- 
cal particle current, local correlation gij+i) and also the first 
two occupancy probabilities Pq and Pi, are presented as a 
function of the ramp time r. This figure clearly uncovers the 
existence of two dynamical behaviors. First, we observe that 
for short ramp times, the densities at ^ = 32 which lies in the 
center of the trap and / — 18 which lies close to the forming 
Mott-insulating barrier are both approximately constant, fol- 
lowing the evolution of the homogeneous system |69]. In fact, 
variations (and oscillations) of both central and outer densities 
become significant only for longer ramp times, beyond h/ J. 
In contrast, the evolution of both the occupancy probabilities 
and the compressibility occurs on a much faster time-scale: 
these observables vary rapidly at short t and display less pro- 
nounced variations at larger r. These two distinct behaviors 
reveal the presence of two dynamical regimes |17, 39, 40]: 
(i) the intrinsic dynamics, here occurring at short-times be- 
fore the particle transport sets in (present in both the homoge- 
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neous and trapped systems); (ii) a long-time behavior associ- 
ated with particle transport and clearly due to the inhomoge- 
neous structure of the density profiles. 

Qualitatively, one can understand the origin of different 
time-scales from the continuity equations of Sec. Ill CI For 
instance, the incoming and outgoing particle currents balance 
each other in a translational invariant configuration so that the 
density remains constant. However, for the density fluctua- 
tions whose current operator O has correlated terms, no such 
balance is achieved and consequently these quantities evolve 
with time. In the case of an inhomogeneous gas, gradients of 
local quantities and chemical potentials inevitably give rise to 
particle currents, themselves sustaining the evolution of all lo- 
cal quantities. Contrary to the intrinsic dynamics, we expect 
these effects to vanish with reducing the trap amplitude Vq. 
Thus, their time-scale is distinct from the intrinsic one and is 
related to the external potential strength. In order to quantify 
better these ideas, we now present arguments based on pertur- 
bative calculations. 



2. Insights from perturbative expansions 

In RefT?, we observed that, for the quench parameters typ- 
ically considered, the homogeneous dynamics of most local 
observables was well reproduced by time-dependent pertur- 
bation theory, particularly in the small-r regime. Working in 
the initial Hamiltonian eigenstates basis \a), of energy Ea, 
the first-order expansion in 5U /r for a real, symmetric and 
dimensionless observable O reads: 



0(r, 5U) = Oo 



OoalaO 



(21) 

The frequencies Ua = {Ea — Eo)/h are excitation ener- 
gies of level I a) with respect to the ground-state |0). la/s = 
^;(a|/j|/3) are the matrix elements of the interaction opera- 
tor and Oai3 = {a\0\/3) those of the observable. Such series 
can be well-behaved in the thermodynamic limit even in the 
absence of a spectral gap, and this is what we observe for our 
setup by looking at different system sizes. Taking the r 
limit Ii70i1 . the response of the observable is typically quadratic 



0(r, SU) ~ Oo 



(22) 



where the ± sign depends on the observable. We have here 
introduce the "curvature" 



•'^ 3 J Ooo 



(23) 



containing an intrinsic characteristic ramp time tq associated 
with the observable O ( J is there for dimensionality normal- 
ization): 



Tq'^ — — I LUgOaglal 



Q#0 



(24) 



The curvature helps understand the departure from the initial 
value Ooo, one can see, for example, in Fig.[T] In partic- 
ular, fo is linear with the quench amplitude 6U (within this 
approximation). For example, in the homogeneous gas limit, 
one can obtain an explicit expression for the driving of the par- 
ticle fluctuations /„2 by, in addition, resorting to perturbation 
theory in J /Ui (strong interaction limit). We find that 



fn 



32 5UJ^ 



(25) 



which is consistent with the change of the compressibility at 
short time plotted in Fig. [T] The breakdown of the quadratic 
behavior, which coincides with the onset of the relaxation, is 
expected to happen on a time scale t = |/„2 We note 
that the parameter J/Ui in Fig. [T] is not in the regime where 
perturbation theory is expected to give a quantitative descrip- 
tion. Nevertheless, putting numerical values in ( |25l l, one finds 
short relaxation times (below h/ J) compatible with Fig.[T] 

It is also worth mentioning that in the definition of the 
curvature ( |23] ), we were careful to separate what depends 
on the quench protocol, the parameter SU and the prefactor 
2/3, from what is intrinsic to the initial ground-state: Ooo 
and To- Indeed, when the U{t) function is of the general 
type SU /{t/r), the prefactor 2/3 is replaced by 4 dx{l — 
x)f{x). Hence, tq is an intrinsic characteristic time of the 
initial state. We stress that the quantities in ( |23] ) and ( l24l l are 
accessible by ground state numerical techniques. Within this 
perturbative framework, one can easily understand the two- 
regimes discussed above and also derive, in the homogeneous 
case, relations between the characteristic time-scales of vari- 
ous observables. 

We first consider the characteristic time associated with the 
local density operator. In the homogeneous case, the ground 
state is characterized by a spatially uniform local density. Tak- 
ing advantage of this symmetry, we find that r„ = oo. This 
result agrees with the fact that the density remains constant 
for all times. In contrast, for the local density fluctuations (or 
compressibility) and local kinetic energy, r„2 /g are finite even 
in a homogeneous system since the matrix elements in (l24l) do 
not vanish. Furthermore, the two time-scales are actually re- 
lated to each other. Since the Hamiltonian has only two terms, 
we find that ^,(0|gM+i|a) = {U,/J) E/(0|«?|a) leading to 
Tg = ^J^Ty-fi . This relation agrees with a dimensional anal- 
ysis of ([Tol l and is also in qualitative agreement with Fig. [T] 
where we find a slightly slower relaxation for the kinetic term 
as compared to the compressibility. In addition, these time- 
scales are themselves related to the characteristic ramp time 
for the heat, Tc, defined as = y^f-j^ X]q \IaO 
Then, we find that t„2 = Tc/ -\/24 (although the prefactor de- 
pends on the chosen definition for Tc). 

We now turn to the situation where a trapping confinement 
is present. In this case, the translational symmetry is lost 
leading to a finite t„. Naturally, Tr^/g should also be af- 
fected by the presence of the trap, but provided the latter is 
small enough, the corrections can be negligible as illustrated 
by Fig.[T] We expect that T„(Vb) diverges when the trap mag- 
nitude Vq reaches zero. Consequently, by tuning Vo to a low 
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FIG. 2. (color online). Final local density and compressibility pro- 
files after a slow quench from Ui — 4J to Uf = 6J for different 
ramp times, r, and for the ground state (r — oo) at f/ = 6 J in the 
trapped system. Left panels TV = 24, right panels = 48. 
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enough value, one should in general be able to observe the 
intrinsic dynamics of the system occurring below t„. The be- 
havior of the T„(Vb) function is an open issue, particularly 
because Vo is not a perturbation in experiments and realistic 
numerical calculations. If we were to trust a naive first order 
perturbation argument for the relatively unphysical situation 
of a gas in a box and perturbed by a small Vq, one would 
expect a linear scaling of the matrix elements, yielding the 
scaling r„ oc l/y^. However, this scaling only serves as 
an illustration of the above statements. Of course, t„ also de- 
pends on Ui I J. Finally, we may argue that when Vq is too 
large, transport phenomena could eventually hide the intrinsic 
dynamics. 

These results put on firmer grounds the existence of two dif- 
ferent dynamical regimes: one deeply connected to inhomo- 
geneities and controlled by Vb, and the intrinsic one present in 
the homogeneous gas and much less sensitive to Vq. 



3. Profiles and "Mott barriers " 

We detail here the spatial evolution of local quantities. 
In Fig. 12] we present the final profiles for the density and 
compressibility as a function of ramp time. At low filling 
(iV = 24), the shape of these final profiles is well understood 
if one resorts to the arguments presented above: we see that 
for short ramp times the density profiles barely evolve while 
the compressibility changes considerably. For longer r, the 
profiles approach smoothly the final ground state configura- 
tion. For the larger filling N — 48, the evolution is more 
involved. A strong reduction of the compressibility occurs lo- 
cally in regions of filling close to unity already for short ramp 



FIG. 3. (color online). Time-evolution of the local compressibility k; 
and current = {ji,i+i) during a slow quench from Ui — 4J to 

Uf = 6 J in a time r — 7h/ J for A'^ = 48 in the trapped system. As 
the "Mott barriers" are formed, the current in their vicinity weakens. 



times T sa h/J while the formation of pronounced Mott- 
insulating "shells" in the density profile only takes place at 
much longer ramp times, after 5h/J. 

To understand better the complex dynamics at play in the 
presence of regions close to filling one, we show in Fig. [3] 
real-time snapshots of the compressibility and particle current 
for the ramp time t = 7 h/J. The connection between these 
two quantities becomes evident at close inspection. One first 
notices that, once again, the compressibility in regions away 
from unit filling evolves quickly while the flow of atoms to- 
wards the system boundaries takes a much longer time to set 
in. In addition, once the compressibility is sufficiently sup- 
pressed in the regions of filling one, the current in these re- 
gions weakens, which slows down the density redistribution 
across the gas. Even though the regions close to unit filling 
are small and are not real Mott-insulating plateaus, they still 
reduce significantly the transport from the inner to the outer 
superfluid domains. Consequently, the onset of low compress- 
ibility regions explains why systems above unity filling evolve 
slowly when U is increased. From here on, we will refer to 
these regions as "Mott barriers". 

To shed even more light on the build up and suppression of 
the particle current, we analyze the contribution of the differ- 
ent terms appearing in ( fT2] i which make up the time-derivative 
of the particle current. For each term, we plot in Fig. |4] the 
different contributions to — div( J^) for various times in order 
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FIG. 4. (a — d) Time-evolution of the different contributions to the 
time-derivative of the particle current: Eq. J12at (a), Eq. l ll2bb (fe), 
Eq. J12cb (c) and Eq. J12dt (d). (e) time-evolution of the time- 
derivative of the particle current. (/) time evolution of the particle 
current. For (a) to (e), the value at f = is subtracted. The evolution 
parameters are the same as in Fig. [3] 



to understand what drives the evolution of the particle cur- 
rent. The first remarkable feature is Fig.|4|e) where the time- 
derivative of the current changes sign near unity filling around 
t ~ 5h/ J. This inversion is a clear indication that the current 
is being suppressed by the formation of Mott barriers. By con- 
sidering each contribution separately using Fig.HJa-d), we ob- 
serve that the main driving terms boosting the particle currents 
are the ones related to the local kinetic energy ( I12al ) and ( I12cl ). 
and the density assisted hoppings term ( I12db while the density 
gradients ( |12b| ) become significant only at the edges where 
the density varies rapidly. The most striking phenomenon is 
due to the density assisted hoppings term ( I12db . We see on 
Fig. Htd) that this term, which is non-zero only in the pres- 
ence of interactions, changes sign in the regions where Mott 
barriers are forming thus drastically slowing down the equili- 
brating out-flow of atoms. 



We finally stress again that the time-scales associated with 
the contributions ( fT2] l are essentially controlled by the steep- 
ness of profiles induced by the trapping potential Vq. Within 
our choice of parameters they take longer times than the in- 
trinsic evolution. Experimentally, changing the confinement 
strength Vq/ J would affect both in time and magnitude the 
creation of Mott barriers Ii7 1|1 . 
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FIG. 5. (color online). Final local energy (a, e) and local excess en- 
ergy (b, /) profiles after a slow quench from Ui = 4J to Uf — 6J 
for different ramp times, r, and for the ground state (r — oo) at 
U = 6J for a trapped system, (c) is the contribution to qi due to 
the external operator (see ( |19cb ) while (d) is the contribution from 
the energy currents (see l |19bb ). The dashed line in (6) and (/) cor- 
responds to {hi)o.i — {hi)o.f- 



4. Energy transport and heat production 

We now turn our attention to the energy transport and heat 
production during a quench. We present in Fig.|5ja, e) the fi- 
nal local energy profiles for N — 24 and 48 |72]. This figure 
confirms our findings obtained from the analysis of the local 
density and compressibility profiles: for = 24 the system 
approaches the adiabatic limit much faster than for N = 48. 
For N = 48, the final profile remains highly excited even 
for the longest ramp time considered (r « 25ft/ J). The local 
excess energy production highlights a series of differences be- 
tween the evolution of systems with filling below and above 
one (see Fig. |3b) and Fig. |3f)). We first notice that the lo- 
cal excess energy is smaller by nearly an order of magnitude 
for iV = 24 compared to = 48. Furthermore, while for 
r = 25ft/ J and = 24 the local excess energy is rather 
uniformly distributed and close to zero, the iV = 48 result ex- 
hibits strong spatial fluctuations with qi large and negative at 
the edges and large and positive at the center of the cloud. In 
fact, the local excess energy pattern resulting from the quench 
is highly non-trivial even for the seemingly simplest situation 



9 



-0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 




1 30 60 

site / 



FIG. 6. Time-evolution of the different contributions to the energy 
current: (a) Eq. (Ha), (b) Eq. i flibt . (c) Eq. iUc[ and (d) Eq. l fT4dt . 
Time evolution of the full energy current (e). The evolution parame- 
ters are the same as in Fig. [3] 



where TV — 24, as illustrated in Fig.|5jb, c, d). At short ramp 
times (sudden quench limit), particles and energy currents are 
negligible so that the term ( |19bt does not contribute, all the 
final excess energy being a balance between the ground state 
energy difference ( |19a| i and the density fluctuations average 
( |19c| ). The latter is always positive and distributed rather uni- 
formly in a Gaussian-like function whose maximum decreases 
with T (see Fig. |5jc)). Hence, for short ramp times, the bulk 
retains most of the local excess energy while the edges have 
negative qi due to the term ( |19at . For longer ramp times, en- 
ergy currents set in with the effect of redistributing energy 
from the bulk to the edges (see Fig. lUb)). Thus, these cur- 
rents tend to strongly reduce both the spatial fluctuations and 
the total excess energy (heat) produced by the quench. For 
intermediate times, either negative or positive qi at the edges 
and in the bulk (see for instance the opposite distribution for 
r = ion/ J and r 15h/J for N = 24) can be found. This 
effect arises as, for these parameters, the density profiles over- 
shoot their final ground state configurations. For N — 48, the 
"Mott barriers effect" tends to freeze the excess local energy 
pattern to the sudden quench typical distribution with negative 
qi at the edges and positive in the bulk. Let us note that the 
freezing of the local excess energy pattern is strongly related 
to the frozen density pattern. 

Looking at the contributions to the evolution of the energy 
current in Fig.|6] we can identify the leading contribution driv- 
ing the energy redistribution. Comparing Fig.Uf) andjSfe), 
we see that the overall evolution of the particle and energy 
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FIG. 7. (color online). Total heat Q(r) vs. inverse ramp time r for 
both directions of the quenches for N = 24: and A'' = 48. 

current is very similar: both currents set in at about the same 
time, and are suppressed in regions where Mott barriers form. 
We also observe that the main contribution, determining the 
sign, to the energy current is from the density-assisted parti- 
cle current (see Fig. |6fd)). However, this flow of energy to- 
wards the system edges is partially counterbalanced by two 
terms (Fig.|6|b) and (c)) where the energy transport occurs in 
the opposite direction to the particle current. 

Finally, one may wonder how the total heat produced Q as a 
function of 1 /r differs from the homogeneous situation stud- 
ied in Ref. '17'. We show how the heat behaves as a function 
of the ramp time in Fig. |7]for iV = 24 and 48 for both the 
quench from Ui = 4 J to C// = 6 J and its reverse. We first 
notice that these curves are more complex than the one pre- 
sented in Ref.[l7|for a homogeneous system. We also observe 
that the heat per atom produced in the case = 24 is always 
much lower than the one for N — 48. Our understanding of 
this phenomenon is that for lower filling the populated excited 
states are less energetic as they are less likely to have doubly 
occupied sites. We finally observe that, for fast ramps, the heat 
produced in the protocol with Ui = 4 J and Uf = 6 J is larger 
than for the reverse protocol (while the opposite happens for 
slower ramps). We relate this to the fact that in the Ui/J — 4 
initial state a lot of particle fluctuations are present leading to 
a large interaction energy in the final state. However, in or- 
der to fully understand the crossover to the inverse behavior 
at slow ramp times, the number of excitations that are created 
and their final energies would need to be identified, a task that 
we leave to future studies. 



5. Comparison with mean-field Gutzwiller method 

Our aim here is to understand to what extent the mean-field 
Gutzwiller method can describe the time-evolution of a Bose 
gas loaded into a one-dimensional optical lattice and confined 
to a parabolic trap. With this objective in mind, we study 
here a system made of 54 atoms confined to a parabolic trap 
with Vb = 0.006 J and loaded in an optical lattice of 84 sites. 
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FIG. 8. Time-evolution, within the Gutzwiller method, of a trapped 
one-dimensional Bose gas loaded into an optical lattice during a 
quench from Ui = 6Jtof// — 15 J for two different ramp times. 
Plotted quantities are the local density (fii), the local compressibil- 
ity Hi, the current and the superfluid order parameter (fej). 
TV = 54, L = 84, Vt = 0.006J. Upper panels: r = 16h/J. 
Lower panels: r = 60h/J. 



and consider slow quenches from Ui ~ 6, J to Uf ~ 15 J 
for two different ramp times: r = 60ft/ J and r = 16h/J. 
These quenches begin on the superfluid side and the inter- 
action strength is linearly increased up to a value above the 
71 = 1 homogeneous superfluid-Mott-insulating transition, 
occurring at Uc ~ 11.7 J (using the Gutzwiller method). At 
mean-field level, the ground state at U — 6 J is a superfluid 
with a central density above one while the ground state at 
[/ = 15 J presents a broad Mott plateau. 



Considering Fig. |8] we first notice that the Gutzwiller 
method captures well the presence of two dynamical regimes. 
For both ramps, we see that the evolution of the local com- 
pressibility and superfluid order parameter begins at t ~ 
whereas the local density and the particle cuiTent remain fixed 
to their initial values for a few h/J. For a sufficiently fast 
quench, as shown in the upper panels of Fig. [8] we see that the 
superfluid order parameter, the compressibility and the cur- 
rent are strongly suppressed in a narrow region around filling 
one. The local suppression of these three quantities around 
t = 6H/ J signals the formation of Mott barriers hindering the 
flow of atoms. For r — 16 J/fi, these barriers are unstable and 
we notice the presence of oscillations reminiscent of the ones 
arising when a strongly interacting phase is abruptly quenched 
to strong interactions [|73ll . 

By comparison, for sufficiently slow ramps, a stable Mott- 
insulating plateau forms at long times. On this plateau, the 
condensate order parameter and the compressibility drop to 
zero. This total suppression of the density fluctuations is an ar- 
tifact of the mean-field method and also results in the absence 
of particle cuiTent on the plateau as, within the Gutzwiller pic- 
ture, the current factorizes into ji-ij — 2J Sy((fe;_i)*(6;)). 
Finally, we also observe in the lower panels of Fig. [8] that 
the quench triggers collective breathing modes signaled by 
density oscillations (along the time axis) in boundary re- 
gions 1 74]. 

From this discussion of slow superfluid-Mott-insulating 
quenches within the Gutzwiller method, we conclude that this 
approach captures some of the important out-of-equilibrium 
physical phenomena uncovered by t-DMRG, however as ex- 
pected it cannot provide an accurate quantitative picture. 



B. Evolution of non-local quantities 

1. Real-space correlations 

Local and non-local correlations can propagate very dif- 
ferently during a quench. To understand how correlations 
evolve during the slow quench of a global parameter, we in- 
vestigate here the evolution of single particle coiTelations. 
Past studies on other systems have found that, after a slow 
parameter c hang e, long-range coiTelations take a long time 
to adjust llH I25II . For example, for spin systems described 
by locally acting Hamiltonians, the propagation of correla- 
tions during the slow change of a global parameter was found 
to be bound by a "light-cone". Outside of this light-cone, 
the so-called Lieb-Robinson bound, only exponentially small 
changes to the coiTelations can be detected ITTSllTql . 

We show in Fig. |9] the value of the correlator 
9L/2+i,L/2+i+d at t = T for different ramp times and two 
fillings. The first striking result emerging from our study is 
that the evolution of this correlator is not monotonic with the 
ramp time. We also find that in all cases the short distance 
correlator responds quickly to the increase of the interaction 
strength, and that even for the fastest quenches the final cor- 
relation function differs considerably from the initial ground 
state correlator. Focusing on the left panel of Fig. |9] we see 
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FIG. 9. (color online). Value of the correlator gi,i+d with / = L/2+1 
after a slow quench from Ui — 4 J to Uf = 6 J in a trap for various 
ramp times r and the final ground state (r = oo). Left panel; — 
24. Right panel: N = 48. 
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FIG. 10. (color online). Final value of the interference pattern N{k) 
after a slow quench from f/i = 4 J to t// = 6 J in a trap for different 
ramp times r and for the final ground state (r = oo). Left panel: 
TV = 24. Right panel: iV = 48. 



that at low filling (N = 24) the short distance correlations 
reach their final ground state values for almost all considered 
ramp times. In contrast, the longer range correlations take 
much longer to reach their corresponding ground state values. 
For example, for r — 5h/ J, the long distance correlations 
have clearly not yet relaxed to their final ground state values. 

In the situation where regions with filling above one are 
present (right panel of Fig. |9]l, the evolution is even more 
involved. In this case, the correlator at t — t varies non- 
monotonically with distance and takes negative values for in- 
termediate ramp velocities. Even for the slowest ramps con- 
sidered, the correlator deviates considerably from its final 
ground state value at all distances. Finally, let us note that 
the final ground state correlations present a dip at a distance 
corresponding to the location of regions of filling one. 



2. Interference pattern for experiments 

Part of the complex dynamical behavior presented above 
can be observed experimentally in the time-of-flight interfer- 
ence pattern N{k) defined in (l2]i. As shown on the left panel 
of Fig. [To] at low filling, the interference pattern present a 
peak at fc = which changes in amplitude non-monotonically 
with the ramp time. This behavior reflects the non-monotonic 
variation of long range correlations discussed in the previous 
section. The final interference pattern is very different in the 
presence of regions close to filling one. In this case, a peak at 
fc 7^ develops at intermediate ramp times (see the right panel 
of Fig. [Tol l. This peak signals the strong out-of-equilibrium 
character of the state formed during the slow quench. How- 
ever, the absence of such a peak cannot be used to conclude 
that the system evolves adiabatically. Unfortunately, the inter- 
ference pattern is not as sensitive to out-of-equilibrium fea- 
tures as correlation functions are: N{k) can be dominated 
by large "in-equilibrium" contributions coming from the short 



range correlators. 

IV. MELTING OF MOTT-INSULATING REGIONS 

In this section, we consider a linear quench from Ui = 6 J 
io Uf — 4 J. At Ui, the system presents a sizable Mott- 
insulating "shell" and a superfluid center, while Uj is close to 
the homogeneous superfluid-Mott-insulating transition point. 
The ground state density and compressibility profiles at Uf 
show none of the features associated with the presence of Mott 
regions. Here again we focus on the different aspects of the 
dynamics: time-scales, particle transport, energy production, 
and experimental signatures. 

A. Existence of two dynamical regimes 

When the interaction strength is lowered, the dynamics 
at play are also characterized by "two dynamical regimes". 
However, as seen on Fig.[TT] in this case, the atoms are mov- 
ing towards the center of the system not towards the edges. 

B. Density, compressibility and energy profiles 

Considering the density and compressibility profiles for dif- 
ferent ramp times shown in Fig. [12] we notice that for ramp 
times of the order of 5h/ J the Mott-insulating regions are al- 
most fully melted and that the system is more compressible. 
For example, local dips, initially present, have completely dis- 
appeared and only plateaus remain. For even longer ramp 
times, the final density and compressibility profiles resemble 
closely the U f ground state. Density redistribution occurs at 
a much faster pace when Mott-insulating regions are melted 
away than when they are formed since in the former case Mott 
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FIG. 11. (color online). Slow quench from Ui = 6 J to Uf = 4 J, 
A'^ = 48. Evolution of local observables in the presence of a trap as 
a function of the ramp time t, and compared with that of a homo- 
geneous system (open symbols) having the same initial local density. 
Observables are the density ni, compressibility kj, occupancy proba- 
bilities Po and Pi , neighboring correlation gi,i+i and particle current 



barriers are no longer effective. To conclude this comparison 
between the two protocols, it is interesting to note that, when 
the interaction strength is lowered, energy is transferred from 
the edges to the center of the system as atoms pile up in the 
central region. The opposite occurs when the interaction is 
increased. 



C. Real space correlations and interference patterns 

Even though for long ramp times the density and compress- 
ibility profiles seem to evolve almost adiabatically, the single 
particle correlator 3^/2+1. L/2+i+d indicates that the system 
is still far from equilibrium. As we can see on Fig. [T3la). this 
coiTelator is negative at large distances and remains far from 
its ground state U f value even at long ramp times, except for 
short distances. Therefore, to judge if the system has reached 
equilibrium by solely considering the density and compress- 
ibility profiles is inadequate. Our results show unequivocally 
that the system is far from having fully relaxed even at long 
ramp times. The non-equilibrium nature of the final state can 
be partially probed by measuring experimentally the interfer- 
ence pattern (see dU). On Fig. [13] (right), we see that for in- 
termediate ramp times, the peak at fc = is shifted to higher 
momentum signaling the non-equilibrium nature of the final 
state. However, as the interference pattern is a sum over all- 
distance correlators and is dominated by short-distance val- 
ues, it is difficult to distinguish the shifted peak at long ramp 
times. 



V. CONCLUSION 

In this article we investigated the dynamics of the Mott- 
insulating regions of a bosonic gas trapped and loaded into 
an optical lattice as the interaction strength is changed lin- 
early with time. We considered two situations: we first stud- 
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FIG. 12. (color online). Final profiles for the local density (a), local 
compressibility (6), local energy (c) and local excess energy (d) after 
a slow quench from Ui = 6J to Uf — 4J for different ramp times, 
r, and for the ground state (r — oo) at (7 = 4 J for a trapped system, 
(e) is the contribution to qi due to the external operator (see l |19c| l) 
while (/) is the contribution from the energy currents (see l |19b| (). 
The dashed line in (d) corresponds to {hi)o,i — {hi)oj. 



ied how Mott domains are formed by ramping up the interac- 
tion strength from Ui = 4J to Uf — 6 J and, in a separate 
set of simulations, investigated how the domains melt when 
U is ramped down. We conducted this study by examining 
how the atomic density and compressibility profiles evolve, 
how particles and energy flow through the system, how heat 
is produced and how single particle correlations propagate 
as a function of the ramp time. For both situations we con- 
firmed the existence of two dynamical regimes: an intrinsic 
regime occurring at short times before particle transport sets 
in, and a long time behavior connected to the system inho- 
mogeneities and controlled by the strength of the underlying 
trapping potential. We were able to establish the existence 
of these regimes on firmer grounds using various arguments 
based on time-dependent perturbation theory. In a system with 
regions above unity filling, we found that a linear increase of 
the interaction strength is accompanied by the formation of 
Mott insulating barriers which hinder the flow of atoms from 
the center towards the edges. The emergence of these barriers 
is evidenced by dips in the local compressibility and by the 
suppression of the particle current in regions where the local 
density nears unity. We also established that, in these regions, 
the change in sign of the particle cuiTent time-derivative is due 
to density assisted hopping, a mechanism which only exists 
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FIG. 13. (color online). Ramp from Ui = 6J to Uf = AJ, N ^ 
48. Left panel: final value for correlator gi^i+d with I = L/2 + 1 
for different ramp times, r, and the final ground state (r — oo). 
Right panel: final value for the interference pattern N{k) (see l[2]l) 
for different ramp times and the final ground state (r = oo). 

when U is finite. The presence of theses barriers has mul- 
tiple consequences, among others, the system "freezes" into 
a highly excited configuration and long range single particle 
correlations deviate strongly from their final ground state val- 
ues, even for the slowest ramp considered. This last feature 



could possibly be detected experimentally from the gas inter- 
ference pattern. By comparison, when the interaction strength 
is ramped down the evolution is much less involved. For suf- 
ficiently long ramps, the density, compressibility and local 
energy profiles approach the corresponding Uf ground state 
configuration. However, even for the slowest ramp consid- 
ered, the final system is still far from being equilibrated as the 
associated one-body correlator departs strongly from its final 
ground state value for all distances. To conclude, we believe 
that this thorough investigation of the dynamics of a strongly 
interacting bosonic gas will help experimentalists devise pro- 
tocols to prepare complex quantum phases and provides a new 
perspective to the understanding of the coherent evolution of 
quantum systems in the presence of inhomogeneities. 
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Appendix A: Useful commutators 

We list below some helpful commutators which gener- 
ate the terms appearing in the various current operators of 
SecHTC] 

[buhl] ^ bi 

[blbiM=blbi 
[blbuh^i] = {l + 2hi)b{M 

[bjh,blbi] ^ ^blbk (m^k) 

[blbk, hlbi] ^hi~ hk 

Permutations between / and k are obtained by talcing the her- 
mitian conjugate. 

Appendix B: Equation of evolution for the occupancy 
probability 

A single-site I of the Bose-Hubbard model is fully char- 
acterized by the occupancy probabilities Pm of having ni 
bosons onsite. The reduced density-matrix of the site is diago- 
nal because of the conservation of the total number of bosons 
and reads 

N 

pi{t)=J2PnMni){ni\. (Bl) 

The mean-value and standard deviation of the P„, (t) distri- 
bution are simply {ni){t) and -\/ (k;) (t). In order to get the 
continuity equation for (t), we introduce the characteristic 
function 

p=0 ^' ni 

(B2) 



15 



such that (nf)(i) = {-i)P C| . Using (|B2]i, the proba 



bilities are recovered using 



61=0 



PnM = ^f dOe-^'^'^fie-t). (B3) 
Using the relations 6ie"^"' = e"^("'+i)6; and feje'^"' = 



,z(«,-i)gt^^g get 



(B4) 
(B5) 



By taking the p derivatives of this equation with respect to 
9, we get the time-evolution of (n^). In particular, we recover 
the time-evolution of the local density with the first derivative. 
We also see that this equation yields correlated currents of the 
form lb\bi+i'nFi~^) for the evolution of (nf ). If we want the 
time-evolution of the probabilities (t), we have to integrate 
(IB5l l over 0. Formally, we have: 

dtPn, = ^ r dee'^'^^dtie^''^') . (B6) 

^TT Jo 

While these formulas are of little help for the numerics, they 
highlight the connection between the evolution of the P„, dis- 
tribution and the correlated currents. 

Furthermore, knowing the evolution equation of P„, (t) al- 
lows one to obtain the evolution equation for the associated 



onsite entropy of particle fluctuations 

Sl{t)^-kBY,PnAt) In Pn, (t) (B7) 

This entropy is rigorously defined also in the non-equilibrium 
regime. Indeed, in a superfluid regime, or even for free bosons 
where P„i is Poissonian, many n have significant weights 
leading to large si, while the n = 1 Mott regime is such that 
there is a strong peak at n = 1 with shoulders at n — 0,2, 
associated with a much smaller entropy. Thus, si is sensitive 
to the nature of the local phase/domain. Formally, we imme- 
diately get the equation of evolution from 

dtsi = -fcB^(atP„,)lnP„, . (B8) 



